!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!! This module provides some basic symbolic representations for
!! multivariate polynomials.
!!
!! \n
!!
!! There are at least two possible representations for polynomials: as
!! a collection of monomials in a one-dimensional array or via the
!! coefficients in the \f$d\f$-dimensional Pascal triangle, \f$d\f$
!! being the number of variables. Here, the first option is chosen as
!! the preferred representation, and a conversion is made to the Pascal
!! triangle form whenever it simplifies an operation, for instance the
!! reduction of the representation (see later).
!!
!! For most of the usages of this module, the precise representation
!! is not a concern and the following subset of the module interface
!! should suffice:
!! \li construction of \c t_sympol objects: \c new_sympoly,
!! <tt>assignment(=)</tt>, \c all_sympoly, \c nll_sympoly
!! \li reduction (see later): \c red_sympoly
!! \li algebraic operations: \f$+\f$, \f$-\f$, \f$*\f$ and \f$**\f$
!! (most of these operations are overloaded to work with both
!! <tt>real(wp)</tt> and \c t_sympol operands)
!! \li evaluation: \c ev
!! \li variable transformation: \c var_change
!! \li derivation: \c pderj
!! \li integration: \c me_int and \c re_int
!! \li visualization: \c show
!! \li consistency checks: \c test (and related error codes)
!!
!! In particular, one should not be directly concerned with the fields
!! of the type \c t_sympol.
!!
!! The representation as monomial collection is not unique, because
!! the monomials can be repeated and their order is arbitrary. This
!! fact indeed simplifies many operations, in particular the sum of
!! monomials, which is simply the concatenation of the respective
!! vectors of monomials. However, even rather simple operations on
!! polynomials such as sums, products and exponentiations can lead
!! very quickly to an explosion of the number of monomials, and hence
!! of the computational cost. For this reason, it's better to reduce
!! polynomials very often. In its reduced form, a polynomial contains
!! monomials that can not be summed, and thus contains the minimum
!! number of monomials. Since reducing a polynomial has some overhead,
!! the type \c t_sympol includes a logical field \c reduced indicating
!! whether the polynomial is already reduced. Reducing an already
!! reduced polynomial has almost no computational cost, so that it's
!! safe to reduce it as often as possible.
!!
!! Concerning integration, there two kinds of definite integrals that
!! are often useful in finite element computations: integrals over the
!! unit regular simplex \f$\mathcal{E}_d\f$, which is the regular
!! (i.e. with equal sides) simplex of \f$\mathbb{R}^d\f$ inscribed in
!! the unit sphere, and integrals over the canonical simplex
!! \f$\hat{K}_d\f$, or master element, which is the simplex of
!! \f$\mathbb{R}^d\f$ with vertices
!! \f$(0,\ldots,0),\,\,(1,0,\ldots,0),\ldots,(0,\ldots,0,1)\f$ (see
!! mod_master_el). The first integrals are useful in the computation
!! of symmetric quadrature formulas, because \f$\mathcal{E}_d\f$ is
!! the simplex with the highest symmetry, while the second kind of
!! integrals are used in the definition of finite element bases. Here,
!! the two integrals are implemented in the functions \c re_int and \c
!! me_int, respectively. \bug The function \c re_int still needs to be
!! implemented; the best solution should be a change of variables to
!! the master element followed by a call to \c me_int.
!!
!! We summarize now some implementation details. The representation of
!! a \f$d\f$-variate polynomial of degree \f$k\f$ in the Pascal
!! triangle form is a one-dimensional real array containing the
!! coefficients \f$\alpha_{i_1\ldots i_d}\f$ of the polynomial in
!! reverse lexicographic order
!! \f{displaymath}{
!!   \begin{array}{rcl}
!!   p_k(x_1,\ldots,x_d) =&& \alpha_{0\ldots0}
!!  + \alpha_{10\ldots0}\,x_1 +\ldots+ \alpha_{k0\ldots0}\,x_1^k \\\
!! & +& \alpha_{010\ldots0}\,x_2 + \alpha_{110\ldots0}\,x_1x_2
!!      +\ldots+ \alpha_{k10\ldots0}\,x_1^kx_2 \\\
!! & +& \alpha_{020\ldots0}\,x_2^2 + \alpha_{120\ldots0}\,x_1x_2^2
!!      +\ldots+ \alpha_{k20\ldots0}\,x_1^kx_2^2 \\\
!! & +& \ldots \\\
!! & +& \alpha_{0k0\ldots0}\,x_2^k + \alpha_{1k0\ldots0}\,x_1x_2^k
!!      +\ldots+ \alpha_{kk0\ldots0}\,x_1^kx_2^k \\\
!! & +& \alpha_{0010\ldots0}\,x_3 + \alpha_{1010\ldots0}\,x_1x_3
!!      +\ldots+ \alpha_{k010\ldots0}\,x_1^kx_3 \\\
!! & +& \ldots \\\
!! & +& \alpha_{0k\ldots k}\,x_2^k\ldots x_d^k + \alpha_{1k\ldots k}\,x_1x_2^k\ldots x_3^k
!!      +\ldots+ \alpha_{kk\ldots k}\,x_1^kx_2^k\ldots x_3^k.
!!   \end{array}
!! \f}
!! This representation is rather cumbersome; however there are at
!! least two simple ways of obtaining it in practice:
!! \li storing the coefficients of the polynomial in an array \f$a\f$
!! with \f$d\f$ indexes ranging from 0 o \f$k\f$
!! \f{displaymath}{
!!  a(\underbrace{0:k,\ldots,0:k}_{d\,\,{\rm times}}),
!! \f}
!! such that \f$a(i_1,\ldots,i_d)\f$ is the coefficient of
!! \f$x_1^{i_1}\ldots x_d^{i_d}\f$, the Pascal triangle form is
!! obtained as <tt>reshape( a , (/(k+1)**d/) )</tt>;
!! \li directly placing the coefficient of the monomial
!! \f$x_1^{i_1}\ldots x_d^{i_d}\f$ at position \f$\sum_{j=1}^d
!! i_j(k+1)^{j-1} + 1\f$ (for indices starting from 1).
!!
!! \note This module serves as a unified interface to symbolic
!! polynomial objects, i.e. to the contents of the two modules
!! mod_sympoly and mod_symmon. The module mod_symmon should not be
!! accessed directly, but only through the interface provided by
!! mod_sympoly.
!!
!! \note Array automatic variables of derived type may cause some
!! problems, both because they tend to blow the stack and because they
!! often hit compiler bugs. So, local array variables of derived type
!! are always declared here to be allocatable, even if their dimension
!! is known trivially in terms of the input arguments.
!<----------------------------------------------------------------------
module mod_sympoly

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_symmon, only: &
   mod_symmon_constructor, &
   mod_symmon_destructor,  &
   t_symmon,       &
   symmon,         &
   assignment(=),  &
   operator(+),    &
   operator(-),    &
   operator(*),    &
   operator(**),   &
   ev,             &
   pderj,          &
   show

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_sympoly_constructor, &
   mod_sympoly_destructor,  &
   mod_sympoly_initialized, &
   t_symmon,     &
   t_sympol,     &
   new_sympoly,  &
   nll_sympoly,  &
   all_sympoly,  &
   red_sympoly,  &
   assignment(=),&
   operator(+),  &
   operator(-),  &
   operator(*),  &
   operator(**), &
   ev,           &
   var_change,   &
   pderj,        &
   me_int,       &
   show,         &
   test,         &
   ! error codes
   wrong_input,    &
   wrong_previous, &
   wrong_monomials

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public and private members
 integer, parameter :: &
   wrong_input     = 1, &
   wrong_previous  = 2, &
   wrong_monomials = 3

 !> symbolic polynomial
 !! \note \c nmon can be zero, in which case \c mons is allocated to
 !! size 0 and deg is 0. This is in fact the simplest way to represent
 !! a zero.
 !! \bug \c nmon should be a derived type parameter
 type t_sympol
   ! main fields
   integer :: nmon                        !< # of monomials
   type(t_symmon), allocatable :: mons(:) !< monomials
   integer :: deg                         !< degree
   logical, private :: reduced = .false.
   !> error indicator, used in \c show and \c test
   integer, private :: ierr = 0
 end type t_sympol

 !> Symbolic polynomial in Pascal triangle form.
 !! \note In this representation, an empty \c t_sympol corresponds to
 !! <tt>nsv.eq.0</tt> and <tt>coeffs(1).eq.0</tt>.
 !! \bug \c deg and \c nsv should be derived type parameters
 type t_psympol
   private
   ! main fields
   integer :: degx !< maximum degree for a single variable
   integer :: nsv  !< # of variables
   real(wp), allocatable :: coeffs(:) !< coefficients ( (degx+1)**nsv )
   integer :: ierr = 0
 end type t_psympol

! Module variables

 ! public members
 logical, protected ::               &
   mod_sympoly_initialized = .false.

 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_sympoly'

 interface assignment(=)
   module procedure real2sympoly
 end interface

 interface operator(+)
   module procedure add
 end interface

 interface operator(-)
   module procedure sub, minus
 end interface

 interface operator(*)
   module procedure mult, mult_scal
 end interface

 interface operator(**)
   module procedure pow
 end interface

 interface ev
   module procedure ev_scal, ev_1d
 end interface

 interface pderj
   module procedure pderj_pol, pderj_pol_k
 end interface

 interface me_int
   module procedure me_int_p, me_int_m
 end interface

 interface var_change
   module procedure var_change_pol, var_change_mon
 end interface

 interface show
   module procedure polyshow
 end interface 

 interface test
   module procedure polytest
 end interface 

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_sympoly_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_sympoly_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   else
     call mod_symmon_constructor()
   endif
   !----------------------------------------------

   mod_sympoly_initialized = .true.
 end subroutine mod_sympoly_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_sympoly_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_sympoly_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   call mod_symmon_destructor()

   mod_sympoly_initialized = .false.
 end subroutine mod_sympoly_destructor

!-----------------------------------------------------------------------
 
 !> This is the main function to generate \c t_sympoly objects from
 !! outside this module (other useful functions for generating such
 !! objects are <tt>assignment(=)</tt> and \c all_sympoly). The input
 !! arguments are explained in details in the following; essentially
 !! they describe a collection of monomials. The monomials can be in
 !! any order and there can be repetitions.
 pure function new_sympoly(coefficients,exponents) result(p)
  !> \c coefficients are the coefficients of the monomials.
  real(wp), intent(in) :: coefficients(:)
  !> \c exponents are the exponents of the monomials organized in a
  !! two-dimensional integer array as follows. The number of rows
  !! defines the number of variables of the polynomial, with row
  !! \f$j\f$ giving the exponent of \f$x_j\f$. Set a column to zero to
  !! indicate a constant. The number of columns must match the size of
  !! \c coefficients.
  integer, intent(in) :: exponents(:,:)
  type(t_sympol) :: p

  integer :: i
  type(t_symmon), allocatable :: mons(:)

   if(size(exponents,2).ne.size(coefficients)) then
     p%ierr = wrong_input
     return
   endif

   allocate(mons(size(coefficients)))

   do i=1,size(coefficients)
     mons(i) = symmon(coefficients(i),exponents(:,i))
   enddo

   p = red_sympoly(sympoly(mons))
 
   deallocate(mons)
 end function new_sympoly
 
!-----------------------------------------------------------------------
 
 !> This function returns the dimension of the space
 !! \f$\mathbb{P}_k(\mathbb{R}^d)\f$ of \f$d\f$-variate polynomials of
 !! total degree \f$k\f$. For the derivation of this result, see \c
 !! all_sympoly. \note In this function we deal with polynomials of
 !! maximum degree \f$k\f$, while in the representation used in \c
 !! t_psympol we fix the maximum degree in each variable \c degx.
 pure recursive function nll_sympoly(d,k) result(n)
  integer, intent(in) :: d, k
  integer :: n

  integer :: i

   if(d.eq.0) then
     n = 1
   elseif(d.eq.1) then
     n = k+1
   else
     n = 0
     do i=0,k
       n = n + nll_sympoly(d-1,i)
     enddo
   endif
 
 end function nll_sympoly
 
!-----------------------------------------------------------------------

 !> This function returns a basis of the space
 !! \f$\mathbb{P}_k(\mathbb{R}^d)\f$ of \f$d\f$-variate polynomials of
 !! total degree \f$k\f$ composed of single monomials. The cardinality
 !! of the basis is given by <tt>nll_sympoly(d,k)</tt>.
 !!
 !! \n
 !!
 !! The construction of the basis is recursive in the spatial
 !! dimension \f$d\f$ and relies on the fact that, given a basis
 !! \f$\mathcal{B}_k^d\f$ of \f$\mathbb{P}_k(\mathbb{R}^d)\f$, the
 !! higher dimensional basis is computed as
 !! \f{displaymath}{
 !!   \mathcal{B}_k^{d+1} = \left\{ \mathcal{B}_k^d ,\,
 !!     x_{d+1}\mathcal{B}_{k-1}^d,\,
 !!     x_{d+1}^2 \mathcal{B}_{k-2}^d,\ldots,\,
 !!     x_{d+1}^k \right\}.
 !! \f}
 pure recursive function all_sympoly(d,k) result(b)
  integer, intent(in) :: d, k
  type(t_sympol), allocatable :: b(:)

  integer :: i, m, n, expd(d)
  type(t_sympol) :: xd

   ! allocation
   allocate(b(nll_sympoly(d,k)))

   if(d.eq.0) then
     b = 1.0_wp
   else

     ! define the new variable x_d
     expd(1:d-1) = 0; expd(d) = 1
     xd = sympoly((/symmon(1.0_wp,expd)/))

     if(d.eq.1) then
       !b = xd**(/(i,i=0,k)/)
       ! workaround for Fermi
       do i=0,k; b(i+1) = xd**i; enddo
     else
       m = 0
       do i=0,k
         n = nll_sympoly(d-1,k-i)
         b(m+1:m+n) = red_sympoly( (xd**i)*all_sympoly(d-1,k-i) )
         m = m+n
       enddo
     endif
   endif
 
 end function all_sympoly
 
!-----------------------------------------------------------------------
 
 !> Conversion \c t_psympol \f$\mapsto\f$ \c t_sympol. Only terms with
 !! a nonzero coefficients are inserted.
 pure function pascal2sympoly(pp) result(p)
  type(t_psympol), intent(in) :: pp
  type(t_sympol) :: p

  integer :: nmon, i
  type(t_symmon), allocatable :: temp(:)

   allocate(temp(size(pp%coeffs)))

   nmon = 0
   do i=1,size(pp%coeffs)
     if(pp%coeffs(i).ne.0.0_wp) then
       nmon = nmon+1
       temp(nmon) = symmon(pp%coeffs(i),pascal_idx2exp(pp,i))
     endif
   enddo

   p = sympoly(temp(1:nmon))
   p%reduced = .true.
   deallocate(temp)
 end function pascal2sympoly

!-----------------------------------------------------------------------
 
 !> Conversion \c t_sympol \f$\mapsto\f$ \c t_psympol. Since the
 !! Pascal triangle representation is unique, this function can be
 !! used to reduce a polynomial. More precisely:
 !! <tt>pascal2sympoly(sympoly2pascal(p))</tt> is the identity for the
 !! polynomial, but reduces the representation to the minimum number
 !! of monomials.
 pure function sympoly2pascal(p) result(pp)
  type(t_sympol), intent(in) :: p
  type(t_psympol) :: pp

  integer :: i, j
 
   if(p%nmon.eq.0) then ! this representation needs at least one coeff.
     pp%degx = 0
     pp%nsv = 0
     allocate(pp%coeffs(1))
     pp%coeffs(1) = 0.0_wp
   else
     pp%degx = 0
     do i=1,p%nmon
       if(p%mons(i)%nsv.gt.0) &
         pp%degx = max(pp%degx,maxval(p%mons(i)%degs))
     enddo
     pp%nsv = maxval(p%mons%nsv)
     allocate(pp%coeffs( (pp%degx+1)**pp%nsv ))

     pp%coeffs = 0.0_wp
     do i=1,p%nmon
       j = pascal_exp2idx(pp,p%mons(i)%degs)
       pp%coeffs(j) = pp%coeffs(j) + p%mons(i)%coeff
     enddo
   endif
 
 end function sympoly2pascal
 
!-----------------------------------------------------------------------
 
 !> In the Pascal triangle representation, convert a single index into
 !! a vector of exponents (inverse lexicographic order).
 pure function pascal_idx2exp(pp,i) result(idx)
  type(t_psympol), intent(in) :: pp
  integer, intent(in) :: i
  integer :: idx(pp%nsv)

  integer :: n, j, b

   n = i
   do j=pp%nsv,1,-1
     b = (pp%degx+1)**(j-1)
     idx(j) = (n-1)/b ! integer division
     n = n - idx(j)*b
   enddo
 
 end function pascal_idx2exp
 
!-----------------------------------------------------------------------
 
 !> In the Pascal triangle representation, convert a vector of
 !! exponents into a single index (inverse lexicographic order). This
 !! function is the inverse of pascal_idx2exp.
 pure function pascal_exp2idx(pp,idx) result(i)
  type(t_psympol), intent(in) :: pp
  integer, intent(in) :: idx(:)
  integer :: i

  integer :: j, b
 
  ! Careful: sometimes it's useful to pass an idx which size is
  ! smaller than pp%nsv. This happens when working on a polynomial
  ! where some monomial has less variables than the others. So, in the
  ! following loop it's better to use as upper limit size(idx) and not
  ! pp%nsv.
   i = 1
   do j=1,size(idx)
     b = (pp%degx+1)**(j-1)
     i = i + idx(j)*b
   enddo
 
 end function pascal_exp2idx
 
!-----------------------------------------------------------------------
 
 pure function sympoly(mons) result(p)
 ! construct a polynomial from an array of monomials
  type(t_symmon), intent(in) :: mons(:)
  type(t_sympol) :: p
 
   if(any(mons%ierr.ne.0)) then
     p%ierr = wrong_monomials
     return
   endif

   p%nmon = size(mons)
   allocate(p%mons(p%nmon))
   p%mons = mons
   if(p%nmon.gt.0) then ! careful: maxval needs size > 0
     p%deg = maxval(p%mons%deg)
   else
     p%deg = 0
   endif

 end function sympoly
 
!-----------------------------------------------------------------------

 elemental function add(p1,p2) result(p)
  type(t_sympol), intent(in) :: p1,p2
  type(t_sympol) :: p
 
   if( (p1%ierr.ne.0).or.(p2%ierr.ne.0) ) then
     p%ierr = wrong_previous
     return
   endif

   p = sympoly((/ p1%mons , p2%mons /))

 end function add
 
!-----------------------------------------------------------------------
 
 elemental function sub(p1,p2) result(p)
  type(t_sympol), intent(in) :: p1,p2
  type(t_sympol) :: p
 
  real(wp), parameter :: mone = -1.0_wp

   p = p1 + mone*p2
 end function sub
 
!-----------------------------------------------------------------------

 elemental function mult(p1,p2) result(p)
  type(t_sympol), intent(in) :: p1,p2
  type(t_sympol) :: p
 
  integer :: i, j, h
  type(t_symmon), allocatable :: mons(:)

   if( (p1%ierr.ne.0).or.(p2%ierr.ne.0) ) then
     p%ierr = wrong_previous
     return
   endif

   allocate(mons(p1%nmon*p2%nmon))

   do i=1,p1%nmon
     do j=1,p2%nmon
       h = (i-1)*p2%nmon + j
       mons(h) = p1%mons(i)*p2%mons(j)
     enddo
   enddo
   p = sympoly(mons)

   deallocate(mons)

 end function mult
 
!-----------------------------------------------------------------------
 
 elemental function mult_scal(r,p1) result(p)
  real(wp), intent(in) :: r
  type(t_sympol), intent(in) :: p1
  type(t_sympol) :: p
 
   if(p1%ierr.ne.0) then
     p%ierr = wrong_previous
     return
   endif

   p = sympoly(r*p1%mons)

   p%reduced = p1%reduced

 end function mult_scal
 
!-----------------------------------------------------------------------
 
 elemental subroutine real2sympoly(p,x)
  type(t_sympol), intent(out) :: p
  real(wp), intent(in) :: x

   type(t_symmon) :: m(1)

   m(1) = x
   p = sympoly(m)

   p%reduced = .true.

 end subroutine real2sympoly

!-----------------------------------------------------------------------

 elemental function pow(p,n) result(q)
  type(t_sympol) :: q
  type(t_sympol), intent(in) :: p
  integer, intent(in) :: n
 
  integer :: i

   if(p%ierr.ne.0) then
     q%ierr = wrong_previous
     return
   endif

   q = 1.0_wp
   do i=1,n
     ! the number of monomials can easily become very large: reduce
     q = red_sympoly( q*p )
   enddo
 
 end function pow
 
!-----------------------------------------------------------------------
 
 elemental function minus(p) result(q)
  type(t_sympol) :: q
  type(t_sympol), intent(in) :: p

   q = (-1.0_wp)*p
 
 end function minus
 
!-----------------------------------------------------------------------

 !> Reduce a polynomial to its minimal representation, i.e. the
 !! representation involving the minimum number of monomials.
 !<
 elemental function red_sympoly(p) result(q)
  type(t_sympol), intent(in) :: p
  type(t_sympol) :: q
 
   if(p%reduced) then
     q = p
   else
     q = pascal2sympoly(sympoly2pascal(p))
   endif

 end function red_sympoly
 
!-----------------------------------------------------------------------

 !> Compute \f$\int_{\hat{K}_d}p(x)\,dx\f$ for a \f$d\f$-variate
 !! polynomial \f$p\f$. The integral is computed as the sum of the
 !! integrals of the monomials. For the definition of the canonical
 !! simplex (or master element) \f$\hat{K}_d\f$ see mod_master_el.
 !<
 pure recursive function me_int_p(d,p) result(i)
  integer, intent(in) :: d
  type(t_sympol), intent(in) :: p
  real(wp) :: i
 
  integer :: j
  type(t_sympol) :: pr

   ! The polynomial is reduced because the number of terms can grow
   ! very quickly when the degree is large
   pr = red_sympoly(p)

   i = 0.0_wp
   do j=1,pr%nmon
     i = i + me_int(d,pr%mons(j))
   enddo
 
 end function me_int_p
 
!-----------------------------------------------------------------------
 
 !> Compute \f$\int_{\hat{K}_d}m(x)\,dx\f$ for a \f$d\f$-variate
 !! monomial \f$m\f$. The integral is computed recursively exploiting
 !! \f{displaymath}{
 !!   \int_{\hat{K}_d} x_1^{n_1}\ldots x_d^{n_d}\,dx_1\ldots dx_d = 
 !!   \frac{1}{n_d+1}\int_{\hat{K}_{d-1}} x_1^{n_1}\ldots
 !!   x_{d-1}^{n_{d-1}}(1-x_1-\ldots-x_{d-1})^{n_d+1}\,dx_1\ldots
 !!   dx_{d-1},
 !! \f}
 !! where \f$\hat{K}_d\f$ and \f$\hat{K}_{d-1}\f$ are the canonical
 !! simplexes in dimension \f$d\f$ and \f$d-1\f$, respectively.
 !<
 pure recursive function me_int_m(d,m) result(i)
  integer, intent(in) :: d
  type(t_symmon), intent(in) :: m
  real(wp) :: i

  integer :: j, degs(d)
  integer, allocatable :: e(:)
  type(t_sympol) :: p

   ! We have to be careful because m%nsv=size(m%degs) can be smaller
   ! than d, in which case the missing exponents are understood to be
   ! zero. To simplify this functions, we copy m%degs in degs adding
   ! the missing zeros as required.
   degs = 0
   degs(1:m%nsv) = m%degs

   if(d.eq.0) then
     i = m%coeff
   elseif(d.eq.1) then
     i = m%coeff/real(degs(1)+1,wp)
   else
     ! build the (d-1)-variate integrand
     p = 1.0_wp
     allocate(e(d-1)); e = 0
     do j=1,d-1
       e(j) = 1
       p = p - sympoly( (/ symmon(1.0_wp,e) /) )
       e(j) = 0
     enddo
     deallocate(e)
     p = 1.0_wp/real(degs(d)+1,wp) * red_sympoly(   &
       sympoly( (/symmon(m%coeff,degs(1:d-1))/) ) * &
                   ( p**(degs(d)+1) )           )
     ! compute the new integral
     i = me_int(d-1,p)
   endif

 end function me_int_m
 
!-----------------------------------------------------------------------
 
 pure function ev_scal(p,x) result(z)
  real(wp) :: z
  type(t_sympol), intent(in) :: p
  real(wp), intent(in) :: x(:)
 
  integer :: i

   z = 0.0_wp
   do i=1,p%nmon
     z = z + ev(p%mons(i),x)
   enddo
 
 end function ev_scal
 
!-----------------------------------------------------------------------

 pure function ev_1d(p,x) result(z)
 ! first index: different variables
 ! second index: different points
  type(t_sympol), intent(in) :: p
  real(wp), intent(in) :: x(:,:)
  real(wp) :: z(size(x,2))
 
  integer :: i

   z = 0.0_wp
   do i=1,p%nmon
     z = z + ev(p%mons(i),x)
   enddo
 
 end function ev_1d
 
!-----------------------------------------------------------------------
 
 elemental function pderj_pol(p,j) result(djp)
  type(t_sympol) :: djp
  type(t_sympol), intent(in) :: p
  integer, intent(in) :: j

   djp = sympoly(pderj(p%mons,j))

 end function pderj_pol
 
!-----------------------------------------------------------------------

 recursive function pderj_pol_k(p,j,k) result(djp)
  type(t_sympol) :: djp
  type(t_sympol), intent(in) :: p
  integer, intent(in) :: j, k
 
   if(k.eq.0) then
     djp = p
   else
     djp = pderj(pderj(p,j),j,k-1)
   endif
 
 end function pderj_pol_k
 
!-----------------------------------------------------------------------
 
 !> Change of variables: the input are a polynomial \f$p(x)\f$ and an
 !! affine transformation specified by a matrix \f$A\f$ and an array
 !! \f$b\f$ such that \f$x=A\,y + b\f$. The result is
 !! \f$q(y)=p(x(y))\f$. The matrix \f$A\f$ doesn't have to be square.
 !! \note The following constraints apply:
 !! \li <tt>size(a,1).eq.size(b)</tt>
 !! \li <tt>size(a,1).ge.p\%mons(:)\%nsv</tt>
 !!
 !! \note This function could be easily generalized to arbitrary
 !! polynomial transformations \f$x=x(y)\f$.
 !!
 !! \warning To make sure to get the correct behavior, add -assume
 !! noold_maxminloc when using the Intel compiler. This is necessary
 !! because the Intel compiler does not respect the standard without
 !! this flag.
 pure function var_change_pol(p,a,b) result(q)
  type(t_sympol), intent(in) :: p
  real(wp), intent(in) :: a(:,:), b(:)
  type(t_sympol) :: q

  integer :: i, j, d, id, eye(size(a,2),size(a,2)), m
  type(t_sympol) :: pr
  type(t_sympol), allocatable :: y(:), x_of_y(:)

   if(p%ierr.ne.0) then
     q%ierr = wrong_previous
     return
   endif
   if(size(a,1).ne.size(b)) then
     q%ierr = wrong_input
     return
   endif

   pr = red_sympoly(p)

   if(any(pr%mons%nsv.gt.size(a,1))) then
     q%ierr = wrong_input
     return
   endif

   m = size(a,2) ! number of variables y

   ! precompute the polynomials  x1=x1(y), x2=x2(y), ...
   ! compiler-bug: the following line requires -assume noold_maxminloc
   ! when using the Intel compiler.
   id = maxloc(pr%mons%nsv,dim=1) ! specify dim=1 to get a scalar
   if(id.eq.0) then ! p is a constant
     q = pr
   else
     allocate(y(size(a,2)))
     ! build the polynomials y1, y2, ...
     eye = 0; do i=1,size(eye,1); eye(i,i) = 1; enddo
     do i=1,m ! each y(i) is composed of one monomial
       y(i) = new_sympoly((/1.0_wp/),eye(:,i:i))
     enddo
     ! build the polynomials x1(y), x2(y), ...
     d = pr%mons(id)%nsv
     allocate(x_of_y(d))
     do i=1,d
       x_of_y(i) = b(i)
       do j=1,m
         x_of_y(i) = x_of_y(i) + a(i,j)*y(j)
       enddo
     enddo

     ! variable substitution
     q = 0.0_wp
     do i=1,pr%nmon
       ! high order terms can lead to a large number of monomials
       q = red_sympoly( q + var_change(pr%mons(i),x_of_y) )
     enddo
     deallocate(y,x_of_y)
   endif

 end function var_change_pol
 
!-----------------------------------------------------------------------
 
 pure function var_change_mon(m,x_of_y) result(q)
 ! Change of variables in one monomial
  type(t_symmon), intent(in) :: m
  type(t_sympol), intent(in) :: x_of_y(:)
  type(t_sympol) :: q

  integer :: i

  q = m%coeff
  do i=1,m%nsv
    q = q*(x_of_y(i)**m%degs(i))
  enddo

 end function var_change_mon
 
!-----------------------------------------------------------------------
 
 !> Display on \c stdout a \c t_sympol object.
 subroutine polyshow(p)
  type(t_sympol), intent(in) :: p
 
  integer :: i

   write(*,*) ''
   select case(p%ierr)
    case(0)

     do i=1,p%nmon
       call show(p%mons(i))
     enddo

    case(wrong_input)
     write(*,*) 'wrong_input'
    case(wrong_previous)
     write(*,*) 'wrong_previous'
    case(wrong_monomials)
     write(*,*) 'wrong_monomials'
   end select
   write(*,*) ''
 
 end subroutine polyshow
 
!-----------------------------------------------------------------------
 
 !> This function tests the consistency of a t_sympol object,
 !! returning
 !! \li 0 if the test is successful
 !! \li a nonzero error code otherwise, which meaning is given by the
 !! parameters \c wrong_input, \c wrong_previous and \c
 !! wrong_monomials.
 !<
 elemental function polytest(p) result(ierr)
  type(t_sympol), intent(in) :: p
  integer :: ierr

   ierr = p%ierr
 
 end function polytest
 
!-----------------------------------------------------------------------

 !> This subroutine tests the main functionalities of this module. It
 !! should not be used in real application, but it can be useful to
 !! check changes to the module. To use it, add it to the \c PUBLIC
 !! list and simply call it in a program like the following \code
 !! program testmod
 !!  use mod_messages
 !!  use mod_kinds
 !!  use mod_sympoly
 !!   call mod_messages_constructor()
 !!   call mod_kinds_constructor()
 !!   call mod_sympoly_constructor()
 !!   call test_module()
 !!   call mod_sympoly_destructor()
 !!   call mod_kinds_destructor()
 !!   call mod_messages_destructor()
 !! end program testmod
 !! \endcode
 !<
 subroutine test_module()

  integer :: i,j,k
  type(t_sympol) :: f, p(4,4,4)
  type(t_sympol), allocatable :: base(:)

   write(*,*) "Testing the construction of t_sympol objects"
   do i=1,4
     do j=1,4
       do k=1,4
         p(i,j,k) = new_sympoly((/1.5_wp/),reshape((/i-1,j-1,k-1/),(/3,1/)))
         write(*,*) "  monomial of degrees",i-1,j-1,k-1
         call polyshow(p(i,j,k))
       enddo
     enddo
   enddo

   write(*,*) "  monomial basis for d=4, k=5"
   allocate(base(nll_sympoly(4,5)))
   base = all_sympoly(4,5)
   do i=1,size(base)
     call polyshow(base(i))
   enddo

   write(*,*) "Testing integration"
   p = (1.0_wp/1.5_wp)*p
   f = 5.0_wp
   f = ( p(3,1,1)**2-p(1,3,1)**2+p(1,1,3)**2-p(2,3,4) )**2 + p(4,1,1) &
       - f
   write(*,*) "  int( (x^2-y^2+z^2-x*y^2*z^3)^2 + x^3 - 5 ) =", &
     me_int( 3 , f )

 end subroutine test_module

!-----------------------------------------------------------------------

end module mod_sympoly

